#version 430
#define BASIS_OUT 8
#include "eustis_builtin_uniforms"

uniform vec2 theta_phi_range;
uniform float radius;
uniform int layer_id;
uniform int msi_w;
uniform int msi_h;
uniform int msi_layers;
uniform int msi_sublayers;
uniform sampler3D msi_rgb;
uniform sampler3D msi_alpha;
uniform sampler3D basis_k[BASIS_OUT];
uniform samplerCube basis_cubemap[BASIS_OUT];
in vec3 pos_near;
in vec3 pos_far;
out vec4 color;

const float half_PI = 3.1415926535897932 / 2.0;

float ray_sphere_intersection(in vec3 ray_o, in vec3 ray_d, float radius) {
  // (ray_o + ray_d * t - sphere_center) ^ 2 == radius ^ 2
  // At^2 + Bt + C = 0, solve t

  float A = dot(ray_d, ray_d);
  float B = 2.0 * dot(ray_d, ray_o);
  float C = dot(ray_o, ray_o) - radius * radius;

  float det = B * B - 4.0 * A * C;
  if (det < 0.001) {
    discard;
  }
  float sqrt_det = sqrt(det);
  return (-B + sqrt_det) / (A * 2.0);
}

void main() {
  vec3 ray_o = pos_near;
  vec3 ray_d = normalize(pos_far - pos_near);
  float t = ray_sphere_intersection(ray_o, ray_d, radius);

  vec3 p = ray_o + t * ray_d;
  float theta = atan(p.z, p.x);
  theta -= (theta > half_PI) ? half_PI * 4 : 0;
  float n_theta = (theta + half_PI) / theta_phi_range.x + 0.5; // [0, 1]
  float u = n_theta - (n_theta - 0.5) / msi_w;
  float phi = atan(length(p.xz), p.y);
  float n_phi = (phi - half_PI) / theta_phi_range.y + 0.5; // [0, 1]
  float v = n_phi - (n_phi - 0.5) / msi_h;

  vec3 c = vec3(0, 0, 0);
  float d1 = (layer_id + 0.5) / (msi_layers * msi_sublayers);
  c += texture(msi_rgb, vec3(u, v, d1)).xyz;
  float a = texture(msi_alpha, vec3(u, v, d1)).x;

  vec3 ray_d_cubemap = vec3(ray_d.x, ray_d.y, -ray_d.z);
  int msi_layer_id = layer_id / msi_sublayers;
  float d2 = (msi_layer_id + 0.5) / msi_layers;
  for (int i = 0; i < BASIS_OUT; ++i) {
    c += texture(basis_k[i], vec3(u, v, d2)).xyz *
      texture(basis_cubemap[i], ray_d_cubemap).x;
  }
  c = clamp(c, 0.0, 1.0);

  color = vec4(c, a);
}
